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Abstract 

In this study we present evidence that anthropogenic stressors can reduce the resilience of age-structured populations. 
Enhancement of disturbance in a model-based Daphnia population lead to a repression of chaotic population dynamics at 
the same time Increasing the degree of synchrony between the population's age classes. Based on the theory of chaos- 
mediated survival an Increased risk of extinction was revealed for this population exposed to high concentrations of a 
chemical stressor. The Lyapunov coefficient was supposed to be a useful indicator to detect disturbance thresholds leading 
to alterations in population dynamics. One possible explanation could be a discrete change in attractor orientation due to 
external disturbance. The statistical analysis of Lyapunov coefficient distribution is proposed as a methodology to test for 
significant non-linear effects of general disturbance on populations. Although many new questions arose, this study forms a 
theoretical basis for a dynamical definition of population recovery. 
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introduction 

Recovery, resilience and disturbance 

Ecosystems are increasingly affected by human activity. 
Structural as well as chemical impacts threaten ecological 
populations, communities, the connected ecosystem functions, 
processes and services [1]. The questions how the constituents of 
ecosystems react to those impairments and whether they can 
recover are matter of many ecological research studies (e.g., [2], 
[3], [4], [5], [6], [7], [8]). 

The term recovery is often stressed in ecology and monitoring. 
In its broadest sense recovery is defined as the recurrence of a 
system to a preceding status after disturbance, often quantified as 
the time to achieve this pre-disturbance status. Thereby the 
concept of recovery is used to quantify how severely a system is 
affected by a disturbance but also to assess whether a system is 
enduringly negatively affected. 

Recovery constitutes a part of the much debated concept of 
ecological stability which is connected to various other ecological 
aspects as resilience [9], [10]. Like for recovery different 
definitions of the term resilience exist (e.g., [1 1], [12]). Engineering 
resilience is defined as the time that a system takes to return to 
equUibrium following perturbation [13], whereas ecological 
resilience is defined as the amount of perturbation a system can 
withstand before it moves into a different basin of attraction or 
stability domain [14]. Often resilience is simply defined as the 
capacity of the system to return to a given state after a disturbance 
[15]. In this last case recovery and resilience are used synony- 
mously (resilience = rate of recovery). 



As resilience reflects the ability of a community or a population 
to withstand and recover from disturbance it is assumed a 
fundamental property of ecological systems and thus has been 
proposed as an indicator of ecological health and integrity [16]. 
Within this context disturbance is broadly defined as any relatively 
discrete event that disrupts ecosystem, community, or population 
structure and changes resources, substrate availability or the 
physical environment [17]. 

It is generally accepted that most species communities are 
constandy exposed to natural disturbance resulting in changes of 
community structure and life history characteristics of the 
constituting species. From the idea that many of the environmental 
parameters that determine resilience of populations and commu- 
nities to natural disturbance may also influence the response to 
chemical stressors [16] it can be concluded that exposure to 
contaminants or other anthropogenic impacts alters the ability to 
withstand and recover from natural disturbance via reduced 
resilience. 

Two different types of disturbance can be distinguished. On the 
one hand pulse disturbance [18] results in prompt change in 
species abundance within a community (e.g. forest fire or chemical 
spill that reduces densities of certain species). Effects of pulse 
disturbances are generally assumed to recede after some time so 
that the system can return to its initial equilibrium state. On the 
other hand press disturbance causes permanent change in species 
abundance and often leads to the loss of some taxa and 
establishment of an alternative community (e.g. increased 
temperature associated with global climate change or continuous 
input of toxic material into a system). Until the elimination of the 
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stressor a community under continuous press disturbance is 
generally not accepted to return to its original condition [16]. 

Dynamics of age-structured populations 

The fact that age structure is an important feature of naturad 
populations has been deduced in many studies based on 
population models [19], [20], [21]. Many animal species, like 
crustaceans, have life-histories composed of a sequence of stages 
within which their individuals' characteristics are broadly similar 
to those of other individuals in the same stage and markedly 
different from those of individuals in other stages [22]. Within 
these complex life-cycles the vital rates (rates of survival, 
development, and reproduction) almost always depend on 
individuals' age, size, or development stage. It has become widely 
recognized that this diversity can have a significant influence upon 
dynamics at the population level and that in order to understand 
the dynamics of most populations it is necessary to take this 
diversity adequately into account [23], [24]. 

Such variation can be described with the help of time discrete 
nonlinear age-structured population models (for examples see [25], 
[26], [27], [28], [29], [30]) which have been widely studied in the 
biomathematical literature [31] (e.g., [32], [33], [34], [35], [36], 
[37], [38]). A large number of theoretical analyses and numerous 
applications to ecological systems have proven the relevance of 
age-structured population models (e.g. [22], [28], [31], [39], [40], 
[41], [42], [43], [44], [45], [46]). 

Even though these models are of simple structure it was shown 
that they display all kinds of dynamical behaviour, like equilibria, 
periodic cycles, and unpredictable chaotic dynamics [40], [47], 
[48], [49], [50], [51], [52], [53], [54], [55]). The occurrence of tiiis 
complex dynamic behaviour is related to the fact that minor 
changes in parameter or initial values can strikingly change the 
dynamic behaviour of systems [3 1] . Model dynamics were found 
to be very sensitive to changes in life-history parameters as well as 
environmental perturbations [28]. In more general terms of system 
dynamics it was shown that trends and fluctuations in populations 
are determined by complex interactions between extrinsic forcing 
and intrinsic dynamics [46]. 

Daphnia magna as model organism 

Daphnia magna is a Crustacean in the order of Cladocera. In 
terms of assessing the effects of human contaminants D. magna 
exhibits an ideal test organism for this purpose as it constitutes a 
central component in aquatic food web structure [56], [57], [58] 
and is known to be sensiti\'e to a multitude of xenobiotics 
(organophosphates, h(;a\y metals, organochlorines, pyrethroids 
etc.) [59], [60], [61]. Due to its ease of culture and its clonal 
method of reproduction [56], [62] it is commonly used as a model 
organism in ecological risk assessment of pesticides, biocides, and 
industrial chemicals [29]. 

Within a Daphnia population three age classes (roughly equal to 
size classes) representing the different life stages can be 
distinguished, defined by means of age and fertility: neonates (less 
than 24 h old), juveniles (at least 24 hours old but non-mature), 
and adults (mature). These age classes differ in terms of their 
toxicological sensitivity in mortality, reproduction and develop- 
mental time. Effects within the population structure due to 
xenobiotics emerge from these age tiass specific sensitivities, e.g. 
reproduction rate. 

Such a daphnid population exhibits a simple example of an 
ecological system with multi-stage interactions of age-classes a„, 
n = { 1, 2, 3}. Due to its reduced dimensionality and the differential 
interactions within it, the population comes close to the classical 
example of a 3-dimensional dynamical system like the Lorenz 



attractor. This suggests to analyse the effects of disturbance (effect 
= altered intrinsic rate of natural increase) of this system in terms 
of non-linear dynamics. 

DCA as a model substance 

3,4-Dichloroaniline (DCA, CAS-Nr.: 95-76-1, EEC-Nr.: 202- 
448-4) is an intermediate used for the synthesis of plant protection 
products (e.g., linuron, diuron, propanU) and as an azo dye for 
polyester fabrics. Release into the environment mainly occurs 
during use of plant protection agents after biotransformation or 
simply as an impurity of these substances [63]. DCA is known to 
affect survival and reproduction in D. magna already in low 
concentrations, resulting in reduced offspring number and 
increased number of aborted eggs in chronic toxicity tests [64], 
[65], [66]. In higher concentrations DCA has a lethal effect, 
causing increased mortality rates in acute toxicity tests [67]. Thus 
a strong effect of DCA on Daphnia population dynamics is obvious 
[68] and it exhibits a useful model substance to investigate effects 
on population level. 

Non-linear time series analysis 

The autocorrelation function allows for determining the non- 
whiteness of data and detecting periodic components in data. 
Because the autocorrelation function and the power spectrum are 
related by the Fourier transform, they are mathematically 
equivalent, and contain the same information [69]. 

Based on Takens' theorem [70] time-delay embedding allows 
for the reconstruction of the system's phase space and its attractors 
[71], [72], [73]. Scalar observations x(n) and their time delays x(n 
+ kT); k= l,2,...,m— 1 are used to make m-dimensional vectors 
y(n) = [x(n),x(n4-T),...,x(n-l-(m— 1)T) whose components provide a 
coordinate system in which one can identify the attractor structure 
associated with the observations (T = time delay) [74]. 

Recurrence plots are a graphical tool for measuring the time 
constancy of dynamical systems [75]. They represent the 
recurrence of the phase space trajectory to a certain state, which 
is a fundamental property of deterministic dynamical systems [76], 
[77]. The analysis of recurrence plots allows for the detection of 
periodicity and regression [75] as well as typical transitions, e.g., 
bifurcation points occurring in complex systems [78]. 

A Poincare section (also called a surface of section, Poincare 
map or first recurrence map) [79] creates a projection of 
intersection points of an attractor and a hyper-plane in phase 
space [77] . It differs from a recurrence plot in that space, not time, 
determines the plot points. The hyper-plane is adjusted to cut the 
attractor transversal, thereby depicting the orientation of the n- 
dimensional attractor as a cloud of points in (n-l)-dimensional 
space. 

Lyapunov exponents (k) are a measure for exponential 
divergence of initially nearby trajectories in embedded phase 
space and are thus a hint for chaotic dynamics within a dynamical 
system [80]. If a dynamical system exhibits at least one Lyapunov 
exponent larger than zero it is assumed to show sensitivity to initial 
conditions, a prerequisite for chaotic dynamics [81], [82]. 
Calculating maximal global Lyapunov coefficients instead of local 
coefficients has the advantage of averaging out local divergence 
rates [83]. 

Aims 

In this paper the authors present an analysis of the dynamics of 
a model-based age-class structured population of Daphnia magna 
exposed to different concentrations of a chemical stressor (3,4- 
Dichloroaniline, DCA). To support the understanding of the 
effects of chemical disturbance the study focused on the non-linear 
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characteristics of these dynamics and the question whether there is 
evidence that disturbance may lead to non-reversible alterations of 
population dynamics. This was achieved by quantification of 

periodic and chaotic parts of population dynamics, using methods 
from non-linear analysis and concepts from dynamical systems 
theory. 

The hypothesis to be tested in this study was that anthropogenic 
stressors can alter the resilience of an age-structured population 
and thus have a negative influence on the ability to respond to 
further natural disturbance, thereby enhancing the risk of 
population extinction. 

The Lyapunov was tested as an indicator to detect disturbance 
thresholds leading to a discrete change in the investigated 
population dynamics. Interpreting the results in terms of resilience, 
the study was intended to contribute to a dynamical definition of 
population recovery. 

Materials and Methods 

Model-based study 

In princijjk' there are three ways to rc'construct dynamics from 
non-linear systems: (1) analysis of system equations, (2) analysis of 
empirical data and (3) analysis of model-based data. In stress 
ecology simulatitm models are often used to predict recovery time 
[5], [29], [30], [84], [85], [86], [87], [88]. 

Analysis of system equations is more a theoretical task and not 
applicable in this context. Analysis of empirical data is often not 
optimal in the case at hand as the time series available is too short 
to analyze non-linear dynamics sufficiently and detect chaotic 
dynamics. Additionally empirical data are noisy and to find 
deterministic chaos the deterministic skeleton must be extracted 
from the data (which requires long data series as well). Thus, in 
this study we used a validated Daphnia population model (details 
see below) to simulate the empirical data. By this means we were 
able to (1) extend the observational time frame, as the model is 
validated to predict Daphnia population dynamics, (2) extract the 
deterministic skeleton, eliminating the confounding noise (we do 
this by setting chance constant) and (3) take into account the 
diversity within the complex life-cycle to understand the dynamics 
of Daphnia magna [23], [24]. 

Simulated data sets 

The time discrcti; data used in this study were generated by use 
of an individual based population model for Daphnia magna 
(IDamP, [29]). The model predicts the population dynamics 
based on individual life cycles which include feeding rate, growth, 
development, reproduction and survival processes for which food 
conditions and population density (crowding effects) arc the main 
drivers. The life cycle process is represented as descriptive 
regression models, based on a large dataset from D. magna life 
cycle tests. 

IDamP reliably predicts population growth and reproduction 
for D. magna under laboratory conditions using the algae 
Desmodesmus subspicatus as a food source. It was recently shown to 
be able to extrapolate the individual effects of DCA to the 
population level as well as for other toxicants [30], [89], [90]. 
Besides mortality the only sub-lethal effects for DCA considered in 
the model is the inhibition of reproduction. The interactions 
between individuals are due to crowding and competition for food. 
To reduce the model stochasticity all individual properties, except 
a maximum lifetime assigned randomly to the initial individuals 
and the daphnids bom during the simulation, were switched off. 
This maximum lifetime value was not changed anymore during 



the simulation. Raw data from the simulations can be obtained 
from the authors on individual request. 

Five concentration levels for DCA were used in the simulations 
(0 |xg/l, 2 |xg/l, 5 |ig/l, 10 |xg/l, 20 |xg/l, 40 ng/1), held constant 
over the simulation time frame of 365 days thus representing a 
press disturbance. Based on experimental laboratory' data 
concentration-response relationships for DCA on the individual 
level had been calibrated in an integrated toxicity module in the 
IDamP model and validated on population level before [30]. 

As the DCA concentrations are far below the acute mortality 
concentrations for daphnids (48 h LCr,(|: 0.23 mg/1 [63], [67]), it 
can be assumed that the effect wiU mainly occur on the adult 
daphnid reproduction rate. 

A Monte-Carlo set of 100 simulation runs was conducted using 
the same parameter set as published in [30] . Differences in the 100 
simulation runs resulted from differences in initial conditions. 
Before starting the simulations 5 neonate and 3 adult daphnids 
were introduced into the system. These 8 individuals were assigned 
random but fixed values for maximum lifetime. 

As the chaotic parts of the population dynamics would be lost in 
non-linear analysis when based on mean values from these 
simulation runs, non-linear dynamical characteristics were calcu- 
lated for each simulation run separately. Statistics were subse- 
quently calculated from the resulting distributions of non-linear 
characteristics. 

The overall population predicted by the model was additionally 
divided into three age classes (size classes): neonates (<1.4 mm), 
juveniles (>1.4 mm and <2.6 mm) and adults (>2.6 mm). In this 
way the data structure allowed for statistical analysis of population 
age class dynamics. 

Modelling non-linear dynamics 

Data analysis was based on the time series data described above. 
To exclude the effects from transient dynamk:s, the first 50 data 
points were omitted from all analyses. First, 3-dimi;nsional phase 
space was plotted from the raw data of neonate, juvenile and adult 
daphnids. Secondly, exploratory analysis was performed on the 
time series data by use of linear autocorrelation functions. Finally, 
after embedding the time series, non-linear analysis was applied by 
means of recurrence plots, Poincare sectkms and global maximum 
Lyapunov exponents. Data visualisation and analysis was 
performed by use of the software package Tisean [83], [91] in 
MATLAB (Ver. 8.2.0.701, release R2013b) on a linux-based 
system (Debian testing Jessie). The MATLAB source code SI is 
accessible from the supporting information. 

In a 3-dimensional plot neonate, juvenile, and adult Daphnia 
abundances were plotted for each DCA-concentration. Trajecto- 
ries were assessed in terms of compactness and orientation in 
phase space. Autocorrelation analysis was performed with function 
corr from the Tisean package to assess the amount of linear 
structures within the data. To find the best unfolding of the time 
series in the reconstructed phase space a series of time-delay 
embedcUngs was performed by means of the function delay from 
the Tisean package and the optimal parameters were chosen. 
Based on this assessment, finally all subsequent analyses were 
conducted with embedding dimension m = 2 and delay d = 1 
(except the Poincare section, \vliich was performed with m = 3). 
Recurrence analysis [78] was performed to quantify the number 
and duration of recurrences of the system in its state space. We 
used the function recurr from the Tisean package, using m = 2 as 
embedding dimension, d = 1 as delay and = 1 as neighbourhood. 
Poincare sections were created with the function poincare from the 
Tisean package applying an embedding dimension m = 3 and 
delay d = 1 to find a geometric depiction of the trajectories in a 
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lower-dimensional space. Quantification of attractor characteris- 
tics was done by calculation of global Lyapunov exponents. 
Maximum global Lyapunov coefficients for all 1 00 simulation runs 

per age class and concentration level were calculated with the 
function lyap k from the Tisean package [80], [81] using an 
embedding dimension m = 2 and delay d = 1 , revealing a 
distribution of Lyapunov coefficients. Values larger than 1 were 
assumed to be a strong evidence for chaotic dynamics within the 
time series. 

Statistical analysis 

From the distributions of 100 Lyapunov coefficients mean 
values and standard deviations per age class and treatment level 
were calculated and tested for significant differences by means of a 
WUcoxon signed-rank test (function wUcoxon.test in R, two. sided, 
conf level = 0.95, [92]). 

Results 

All results are described for the overall population as well as the 
age classes (neonates, juveniles and adults). Figures are only shown 
for the overall population. The according figures for the age classes 
can be found in the supporting information. 

The mean Daphnia abundance from all 100 simulation runs 
(including the 95% confidence interval) for the overall population 
in the control treatment is exemplarily shown as time series data in 



Fig. 1 (abundances for neonates, juveniles and adults can be found 
in the supporting information. Fig. SI). It was assumed that the 
population and its age classes existed under equilibrium conditions 

under control conditions from t = 50 on. 

Population dynamics were quite comparable for control and 
DCA concentrations lower than 40 |J.g/l (Fig. 2). For the highest 
DCA concentration level the dynamics were obviously difierent. 

In general the autocorrelation functions (Fig. 3) showed that 
with increasing DCA concentration no significant change in 
structure occurred. Comparable periodic dynamics of fluctuating 
correlation coefficients were observed up to 20 |xg/l. This effect 
was especially strong for juveniles and adults, while in the neonate 
age class the effect was clearly less pronounced (see Fig. S2 in 
supporting information). In contrast the autocorrelation function 
for 40 Hg/l showed a very strong correlation structure with clearly 
increased amplitude* and decreased frequency. 

From time-delay embedding (Fig. 4) it was concluded, that the 
trajectories for the overall population settied in the middle of the 
phase space. Especially for n(X)natcs and ju\'eniles (see Fig. S3 in 
supporting information) the trajectories settied within a small part 
of the phase space nearer to the origin for the highest DCA 
concentration (40 Hg/1). 

Increased pattern formation was observed with increasing DCA 
concentrations in the recurrence plots (Fig. 5), especially for 
40 |ig/l. This was also observed for the age classes (see Fig. S4 in 
supporting information). For neonates a very strong lineal pattern 
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Figure 1. Daphnia abundance from all 100 simulation runs for the overall population in control. Solid line: mean abundance, dashed 

line: 95% confidence Interval. 

doi:1 0.1 371/journal.pone.0096270.g001 
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Figure 2. 3-dimensional phase spaces from raw data of simulation run 1 for neonates, juveniles and adults. From top left to bottom 
right, row-wise: control, 2 \ig/\, 5 ng/l, 10 ng/l, 20 ng/l, 40 |i.g/l. 
doi:1 0.1 371/journal.pone.0096270.g002 
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Figure 3. Autocorrelation functions of simulation run 1 for the overall population. From top to bottom: control, 2 (ig/l, 5 (xg/l, 10 (xg/l, 

20 ng/l, 40 ng/l, resp. 

cloi:1 0.1 371/journal.pone.0096270.g003 
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Figure 4. Embedded time series of simulation run 1 for tlie overaii population. From top left to bottom right, row-wise: control, 2 (xg/l, 
5 ng/l, 10 ng/l, 20 ng/l, 40 jig/l, resp. (d = 1, m = 2). 
cloi:1 0.1 371/journal.pone.0096270.g004 



occurred in the highest concentration level, while for juveniles and 
adults this pattern formation was less pronounced. Especially for 
juveniles additional circular structures were observed in the 
highest concentration level. 

Poincare section.s did not show large differences in dynamics for 
treatments from control up to 20 (tg/1 (Fig. 6). For 40 |j,g/l the 
attractor orientation in phase space abruptly changed in location 
(shift towards origin) and cross sectional area decreased. This effect 
was especially pronounced for juveniles and adults (see Fig. 8,5 in 
supporting information). For neonates the attractor was always 
located near the phase space origin. 

Table I shows mean values and standard deviations for 
maximum Lyapunov coefficients calculated from the Lyapunov 
distribution for the overall population in the treatment levels. All 
maximum Lyapunov coefficients for all treatment levels were 
larger than zero (Fig. 7). Lyapunov coefficients of DC A treatments 
up to 20 |Xg/l were only slighdy smaller than for control, while for 
40 ^,g/l a large decrease was observed. 

Table 2 shows the results from the Wilcoxon signed-rank tests 
on statistical differences between mean maximum Lyapunov 
coefficients for all treatment lc\c'ls. The mean maximum 
Lyapunov coefficients for 40 |.lg/l were significantly different from 
all other treatment levels. Also, all mean maximum coefficients for 
the treatment levels were significantly different from the control. 



Discussion 

From the time series data in Fig. 1 it was concluded that under 
control conditions (no disturbance) the population and its age 
classes existed under equilibrium conditions. The existence of 
equilibrium is prerequisite in the attempt to define recovery 
because only if population structure is relatively stable, the 
population can return to defined pre-disturbance conditions if 
given sufficient time following the disturbance [16]. 

Generally the system seemed to converge to a quasi-stable state 
in the long run (see Fig. 2). For concentrations lower than 40 |Xg/l 
the system seemed to setde in a relatively small phase space region. 
This region occurred well bounded (a prerequisite for chaos), but 
the trajectory seemed to follow a complex attractor, exhibiting 
considerable population fluctuations. For the highest DCA 
concentration level fluctuations were much smaller compared to 
lower DCA concentrations. The trajectory setded within a much 
smaller phase space region, nearly approaching the zero abun- 
dance level, especially for the adult age class. These dynamics 
clearly exhibited an increased extinction risk for the overall 
population if this effect would have lasted too long (although this 
does not pose a problem to the age-class itself in the first place, as 
the classes are refilled by aging individuals). 

For the overall population (as well as for all age classes) in the 
control and all treatment levels a sinusoid like oscillation was 
observed (see Fig. 3). Whether unpredictability results from a 
random process or is based on non-linear effects (deterministic 



PLOS ONE I www.plosone.org 



6 



May 2014 I Volume 9 | Issue 5 | e96270 



Chaotic Dynamics in Daphnia Populations 



350 

300 
250 
200 
150 
100 
50 



x\\"* 



350 

300 
250 
200 
150 
100 
50 



» . • . . . " ". • • . 
••- • • 

. •' ' 

»•/ »•* , . * • 

• : • . ./ 
. . ^. : 

:•• 

* •* 

f . 

!*• • 

I 



350 
300 
250 
200 
150 
100 
50 



100 



200 



300 



100 



200 



300 



100 



200 300 



350 
300 
250 
200 
150 
100 
50 
0 



t 



0 



too 



200 



300 



350 
300 
250 
200 
150 
100 
50 
0 



! : i 'Ai'i •!.• ia f 

• ^ . ^ • - - . 



* • V .'V. 



100 



200 



300 



350 




300 








250 




200 








150 




100 




50 




0 





»• -It > 



100 



200 



300 



Figure 5. Recurrence plot of simulation run 1 for the overall population. From top left to bottom right, row-wise: control, 2 |xg/l, 5 |xg/l, 
1 0 ng/l, 20 ng/l, 40 ng/l, resp. (d = 1 , m = 1 ,2, e = 1). 
doi:1 0.1 371/journal.pone.0096270.g005 



chaos) can not be assessed by linear autocorrelation analysis as it is 
unable to distinguish between these two sources of fluctuation. 
However, since individual properties, except random individual 
maximum lifetime, were switched of to reduce the stochasticity of 
the model it is fare assumption that the unpredictability resulted 
exclusively from the non-linear effects. 

For the highest concentration level the amplitude of the periodic 
oscillations clearly increased while the frequency decreased. This 
argues for an increased formation of long-range correlations 
within the data [93] for high concentration levels, enhancing the 
long-term predictability. It was thus concluded that the linear 
predictability of system dynamics increased with increasing 
disturbance. 

As already observed for the raw data phase space, the system 
was forced towards an asymptotically stable equilibrium point for 
all DCA concentrations in embedded phase space (see Fig. 4). The 
fluctuations within the trajectory's settlement area decreased for 
the overall population and especially for the juvenile daphnids for 
the highest concentration level, reflecting the population's 
decreasing ability to respond to changes in environmental 
conditions (decreasing resilience). 

Additionally, as the equilibrium point moved towards the phase 
space origin with increasing DCA concentration, the risk of 
extinction also increased. 



The generally increasing pattern formation in the recurrence 
plots revealed, that the system dynamics increasingly tended to 
return to formerly engaged states when the disturbance increased 
(see Fig. 5). 

While the more uniform spreading of points in the lower 
concentration levels suggested a high portion of chaotic dynamics 
(especially for the neonate daphnids), the lineal pattern found for 
in the highest concentration level suggested a high amount of 
periodic dynamics within the time series [94]. In contrast, the 
circular structures observed for juveniles in the highest concen- 
tration level suggested the occurrence of quasi-periodic dynamics, 
like superimposed harmonic oscillations. These structures can be 
suspected to originate from the realistic and model-built in time 
delay between juvenile response, their maturation and the 
according effects on reproduction. 

For the juvenile daphnids a mixture-like pattern was observed in 
the highest concentration level. It seems like the system settled at a 
transition between quasi-periodic and periodic dynamics. Such a 
transition has already been observed in a population model for 
Triholium [95] , in which the varied parameter of disturbance was 
the rate of larvae developing to adults, resembling the effect on age 
class dependent mortality of juveniles within the Daphnia model. As 
already deduced from the autocorrelation function, from our 
findings it was concluded, that the predictability of the system 
dynamics increased with increasing disturbance. 
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Figure 6. Poincare sections of simulation run 1 for the overall population. From top left to bottom right, row-wise: control, 2 {ig/l, 5 ^g/l, 
1 0 |ig/l, 20 ng/l, 40 ^g/l, resp. (d = 1 , m = 3, C = 1). 
doi:1 0.1 371/journal.pone.0096270.g006 
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Figure 7. Boxplots for maximum Lyapunov coefficients for DCA treatments for the overall population. Boxes show mean value 
(horizontal line), 25 and 75% quantile (box edges), 5 and 95% quantile (whiskers) and outliers>2 sd from mean (circles). 
doi:1 0.1 371/journal.pone.0096270.g007 
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Table 1. Mean values and standard deviations for maximunn Lyapunov coefficients per treatment level for the overall population 
with n= 100. 



Control 2.5 |ig/l 5 (ig/l 10 fig/l 20 )ig/l 40 fig/l 

mean 0.503 0.477 0.470 0.467 0.459 0.261 

stdev 0.0471 0.0492 0.0485 0.0528 0.0507 0.0336 



doi:l 0.1 371/journal.pone.0096270.t001 

The Poincare sections illustrated clearly, that the attractor 
determining the system dynamics was more or less stable from the 
control up to 20 jO-g/l DCA in the overall population as well as all 
age classes (see Figs. 6 and S5). The abruptiy changing location of 
the attractor (esp. for juveniles) suggested that the attractor became 
unstable between 20 |ig/l and 40 Hg/l. 

As all Lyapunov coefficients (5 treatment levels x3 age classes 
xlOO simulation runs = 1500 values) were larger than zero, and 
together with the fact that the dynamics were clearly bounded, 
there was strong evidence for chaotic dynamics within this model- 
based population of Daphnia magna. Consulting the Lyapunov 
coefficient as a measure of rate of chaotic dynamics within the 
population development, this rate decreased for DCA treatments 
compared to control (see Fig. 7). From the biological point of view 
this seems reasonable as positive feedback (influence of exponential 
growth) is likely to be smaller in the presence of a disturbance 
affecting reproduction. 

The decrease of the coeflTicients indicated that higher amounts 
of disturbance decreased the potential for exponential increase of 
trajectory distances. This was already obvious from the former 
analyses, but the statistics on the coeflTicients quantified and 
vahdated this. It suggested some sort of concentration-response 
relationship, thereby approving the Lyapunov coefficient as an 
indicator to detect non-linear effects of population dynamic 
disturbance due to environmental stress. Additionally, the 
significant decrease of the Lyapunov coefficient revealed a true 
and important emergent effect which could not have been 
expected from the raw population abundance pattern. 

Conclusions 

Implications 

The decreasing Lyapunov coefficients with increasing DCA 
concentration levels resembled a concentration-response relation- 
ship. This is an important finding as it is directly connected to the 
important species trait of reproductive capacity that strongly 
triggers the daphnids' population resilience (see discussion below). 
Unfortunately it is not straightforward to quantify the effect as the 



Lyapunov coefficient and the treatment level exhibit a reciprocal 
relationship. Furthermore, the upper effect bound is difficult to 
define, because Lyapunov coefficients approaching zero or values 
below zero mean a discrete change in system dynamics (phase 
transition towards a limit cycle or a stable equilibrium point 
respectively), not a steady change as required for a continuous 
concentration-response relationship. Thus the concentration- 
response relationship would exhibit discontinuities, difficult to 
interpret. In addition, Poincare sections suggested that between 
20 |Xg/l and 40 |ig/l the attractor became unstable and a distinct 
transition to another attractor (basin of attraction) occurred. This 
would mean an extra discontinuity within the system's reaction to 
the increased treatment, which is a further clear contradiction to 
the concept of concentration-response relationships. 

Nevertheless, a no-effect threshold value can easily be inferred 
from the distribution of Lyapunov coefficients. A simple test on 
significant differences between the mean values of controls and 
treatment levels (see Tab. 1) can be conducted. In the present 
study this threshold turned out to be smaller than 2.5 |J,g/l for all 
three age classes as well as the overall population (see Tab. 2). 
Such a threshold could serve as a relevant indicator to assess non- 
linear effects in population time-series data. 

The population dynamics can generally be decomposed into a 
deterministic (periodic plus chaotic) and a stochastic part, the latter 
not being considered in this study. The authors are well aware 
about the diverse discussion about the problem how to distinguish 
stochastic fluctuations from chaotic dynamics (e.g. [96], [97], [98], 
[99], [100], [101]). To reduce sensitivity on stochastic fluctuations 
we used the function lyap k from the Tisean package which is 
said to suppress statistical fluctuations by excluding reference 
points which have no more than nfmin (default = 10) neighbours 
closer than E [83]. 

The control level was found to exhibit a larger amount of 
chaotic dynamics than the treatments. Nevertheless, this does not 
mean that the controls were not predictable and thus not suitable 
for assessment. The dynamics produced by the IDAMP-model 
appeared clearly bounded and represented the naturally occurring 



Table 2. p-values from Wilcoxon signed-rank tests on statistical differences between mean maximum Lyapunov coefficients per 
treatment level for the overall population. 



2.5 [ig/l 5 [ig/l 1 0 \igl\ 20 |ig/l 40 



Control 


4.853E-4 * 


3.489E-6 * 


1.137E-6 * 


2.267E-9 * 


t* 


2.5 jig/i 




0.2401 


0.1183 


6.91 9E-3 * 


f 


5jig/l 






0.5842 


0.09201 


t* 


10 ng/l 








0.316 


t* 


20 ng/l 










t* 



^: p<2.2e — 16, significant differences with p<0.05 are marked by an asterisk. 
doi:l 0.1 371 /journal.pone.0096270.t002 
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fluctuations within an ecologically stable population correctly, esp. 
under undisturbed conditions [30]. 

Natural populations as well as theoretical populations from 
ecosystem models are known to c-xist in a mixture of deterministic 
(periodical, ordered) and chaotic (irregular) dynamics [102], [103]. 
In our study it was shown that the ratio between the two 
deterministic parts of population dynamics (periodic vs. chaotic) 
increased with increasing disturbance resulting in the hypothesis 
that the capability for resilience decreased with increasing 
disturbance. As can be seen from Figs. 4 and 6 the dynamics 
are clearly bounded and the population got more and more 
trapped within a very small phase space region. From these 
findings it can be supposed that the population's flexibility to react 
to further external disturbance was reduced. 

As can be concluded from the Poincare maps used in this study, 
the DCA treatment seemed to destabilize the attractor above 
20 |J,g DCA/1. Consulting the definition of ecological resilience, it 
can be assumed that this effect was caused by altering the phase 
space stability landscape, in the way that the system settled within 
a new basin of attraction. The fact that the system attractor 
became unstable could mean more than the gradually changing 
ratio of chaotic and periodic dynamics with increasing disturbance 
(see autocorrelation function and recurrence plot). It m('ans a 
discrete, qualitative change in system dynamics (phase transition), 
which can be attributed to a shift in demographic parameters (like 
survival or fecundiu ) and has already been reported before for 
insect populations [104]. 

But, generally it is assumed that a press disturbance causes 
sustained alterations within a system's dynamics. A system under 
press disturbance is generally not expected to return to an initial 
configuration [16]. At least for the highest treatment tested it was 
supposed that there was an increased chance, that the system 
would not be able to return to the dynamics typical for less 
disturbed conditions simply by virtue of its own strength when the 
disturbance was removed. Instead it would nc(;d a further 
'disturbance' to force the system back to the attractor which was 
observed for up to 20 ng/1 due to the altered stability landscape. 
Recovery in this context would mean a reconfiguration of the 
stability landscape to a state typical for less disturbed systems. The 
problem about this definition within the experimental context is 
that for treatment levels up to 20 (Xg/l no [:hange of the stability 
landscape was visible at all, whereas for 40 |ig/l a change was 
found indeed. As a reconfiguration step (switching off the press 
disturbance) was not included within the experimented design, it 
was only possible to observe the potential for recovery, not an 
actual reco\'er}' event within the populations. But, from the 
findings in this study it was concluded that the chance for a change 
in attractor orientation increased with increasing treatment level. 

Outlook 

The study revealed that recovery and resilience in age-class 
structured populations exposed to environmental stress are 
complex matters and have to be assessed rather deeply to come 
to relevant conclusions. 

Chaos was observed besides periodicity in the model-based age- 
structured Daphnia population. This was reported before for 
different models, populations under laboratory' conditions and 
even for multi-species communities under laboratory and semi 
field conditions (e.g., [38], [95], [104], [105], [106], [107]). Our 
study connects this fact to the idea of a dynamical definition of 
recovery. In this study the amount of chaotic dynamics turned out 
to be the largest in the undisturbed control treatments, pointing 
out the importance of the natural mechanism of population 
stabilization via fluctuation in terms of resilience, defined as the 



flexibility to changes in environmental conditions (dynamical 
plasticity). The unpredictable part of population dynamics in the 
control treatments showed that the specification of recovery in 
terms of the return to the undisturbed control state maybe a 
questionable definition, as the theory behind the recovery concept 
is usually a static, not a dynamical one. Using Lyapunov 
coefficients as indicators it was shown that a disturbance reduces 
the amount of chaoticity or non-linear dynamics within the 
population. In our opinion this effect gives strong evidence that 
anthropogenic stressors can reduce the resilience of age-structured 
population and foster linear dynamics and resulting periodicity, 
reflecting equilibrium conditions. 

Reduced resilience has a negative influence on populations' 
ability to respond to further disturbances and makes them more 
susceptible to extinction due to climatic changes, for example. This 
is an increasingly important aspect in ecological research which 
has been reported for many systems and their population, some of 
them highly endangered (e.g., [108], [109]). 

These results stand in contrast to some other studies which 
instead suggest the increase of chaos and according decrease of 
stability in populations exposed to environmental stress [110]. But, 
in this study the constant press disturbance in combination with 
the discrete reproduction cycle of the daphnids results in a 
situation comparable to a periodic forcing which has already been 
shown to result in phase transitions in the dynamics of time 
discrete non-linear ecosystem models, either from periodic 
behaviour to chaos [1 1 1] or from chaotic dynamics to periodicity, 
the latter resulting in a chaos control mechanism [112]. 

The study at hand supports the hypothesis that chaotic 
oscillations reduce the degree of synchrony among populations 
or subdivisions of populations (decorrelation effect), thereby 
stabilizing the system and reducing the probability of extinction 
(chaos-mediated sur\'i\al), an interpretation already been pub- 
lished for meta-populations before [113]. Thus, we think that our 
findings add an important aspect to the open discussion about 
resilience and chaotic population dynamics in the presence of 
environmental disturbance, which is stiU an exciting area of 
research [114]. 

The quantitative change in system dynamics and the two 
constituting components can effectively be visualised by linear 
methods (autocorrelation analysis) and non-linear techniques 
(recurrence plot, Poincare section, Lyapunov exponents). Lyapu- 
nov coefficient proved to be a useful indicator to detect 
disturbance thresholds quantifying non-linear effects of contami- 
nants within the Daphnia population. Even a concentration- 
response like relationship was observed for Lyapunov coefficients 
in this study. The modelling of non-linear effects is thus regarded 
as important by the authors, as non-linearities in species 
populations cannot be predicted by linear techniques in the long 
run. 

The authors propose their methodology as a way to test for 
statistically significant non-linear effects of general disturbance on 
dynamics in long-term population data sets, especially when 
assessing the risk that a disturbance (Uke introducing chemicals in 
ecosystems) persistently alters the dynamics. It allows assessing an 
important emergent ecological effect on the population level which 
cannot be retrieved by standard assessment methods, neither 
based on individual endpoints like filtration rates nor based on 
population endpoints like growth rates. 

The aim of this study was not to propose alternative definitions 
for the terms stability or recovery, nor to discuss the appearance of 
non-linear dynamics in population models. To come to a time 
explicit statement about recovery, like in the approach of 
engineering resilience, a larger number of simulation studies have 
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to be conducted, so that a mathematical probability can be 
computed for the fact whether a population dynamic time series 
returns to the control dynamics after a given time. 

Further studies will be done in the future to find a suitable 
statistical representation for a recovery statement, elucidate the 
questions about a lasting altered fitness landscape induced by 
disturbance, calculate mathematical probabilities for this and use 
finer increments in the disturbance concentration levels to detect a 
possible point of attractor instability more precisely. 

Supporting Information 

Figure SI Daphnia abundance from all 100 simulation 
runs. Solid line: mean abundance, dashed line: 95% confidence 
interval. Top panel: neonates, mid panel: juveniles, bottom panel: 
adults. 
(TIFF) 

Figure S2 Autocorrelation function of simulation run 1 . 

From top to bottom: control, 2 |t,g/l, 5 |Xg/l, 10 (ig/l, 20 |ig/l, 
40 Hg/1, resp. Top panel: neonates, mid panel: juveniles, bottom 
panel: adults. 
(TIFF) 

Figure S3 Embedded time series of simulation run 1. 

From top left to bottom right: control, 2 [Xg/l, 5 |Xg/l, 10 |ig/l, 
20 Hg/l, 40 |J.g/l, resp. (d=: 1, m = 2). Top panel: neonates, mid 
panel: juveniles, bottom panel: adults. 
(TIFF) 
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